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We here present the details of the numerical realization of the recently advanced algorithm developed to identify the 
fragmentation in heavy ion reactions. This new algorithm is based on the Simulated Annealing method and is dubbed 
as Simulated Annealing Clusterization Algorithm [SAC A] . We discuss the different parameters used in the Simulated 
Annealing method and present an economical set of the parameters which is based on the extensive analysis carried 
out for the central and peripheral collisions of Au-Au, Nb-Nb and Pb-Pb. These parameters are crucial for the success 
of the algorithm. Our set of optimized parameters gives the same results as the most conservative choice , but is 
very fast. We also discuss the nucleon and fragment exchange processes which are very important for the energy 
minimization and finally present the analysis of the reaction dynamics using the new algorithm. This algorithm is 
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can be applied whenever one wants to identify which of a given number of constituents form bound objects. 
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I. Introduction: 



In recent years, a lot of efforts has been made (in experiments and theory) at low, intermediate and relativistic 
energies to understand the physics which drives heavy ion reactions. A new generation of electronic devices made 
it possible to measure a multitude of observables at the same time which can give information about the hot and 
dense nuclear matter formed during a reaction [jjj . In a heavy ion reaction the density can be as high as 2-4 times the 
normal nuclear matter and one may reach temperatures of about 100 MeV The properties of the nuclear matter 
at high densities are not only of importance for nuclear physics, but are also of great use for the astrophysical studies 
especially for supernova studies. Unfortunately, there is no method to measure directly the properties of hot and 
compressed nuclear matter formed during a reaction Jl[. What one can observe are single hadrons. Their properties 
are mostly determined in the late stage of the expansion and it is quite difficult to find observables sensitive to the 
early stage. For this one has to rely on the theoretical (simulation) models. One can simulate the reaction from the 
start to the end where we find cold nuclear matter in the form of nucleons, light and heavy fragments j^]. The most 
important information which one would like to extract from the simulation are the time scales of different phenomena. 
One would like to know, for example, when particles are created, when fragments are formed, whether they carry any 
information about the hot and compress phase etc.? The key question associated with the time scale of the fragment 
emission is whether it is a thermal or a dynamical process i.e. whether the fragments are created after the system 
has thermalized or can already be recognized early, before a possible thcrmalization sets in . This would point to 
initial-final state correlations Q,?,?,?]. In addition, several conjectures on the equation of state, especially those in 
which the nuclear interaction is strongly momentum dependent, could not be tested so far in simulations because the 
nuclei become unstable. Here an early fragment recognization would allow to study these equations of state. 

We shall concentrate here on multifragmentation. All theoretical models used to study heavy ion collisions are 
based on the transport nucleons and mesons only. Therefore, for the study of multifragmentation a method has to be 
deviced to group the nucleons into free nucleons and fragments. In the past, one has taken the spatial correlations 
among nucleons to group them into fragments ||. Naturally, this approach cannot detect different fragments which 
are (almost) overlapping and therefore will give a single big fragment during the early stage of the reaction where 
density is quite high. In other words, simple coordinate space approaches cannot address the question of the time 
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scale of fragment formation. To study that one needs a method where fragments can be identified even if they are 
overlapping, i.e. methods which are based on phase space. 

We conjecture that in nature at any given moment of the reaction that configuration is realized which gives the 
largest binding energy. That this concept is meaningful and gives sensitive results will be demonstrated later. To find 
the most bound configuration we are confronted with two problems. 

a) The huge number of possible configurations. 

b) The fact that the number of entities changes. Whereas the number of nucleons is constant, the number of free 
nucleons and fragments is a variable. The problem caused by this fact will be discussed later. 

One may approach this problem by simple iterative methods. They, however, do not garranty that a global minimum 
is obtained but may arrive at a local minimum Q . First attempt to overcome this problem has been advanced in ref. 

. Though this method works fine for small systems, its simple numerical implementation poses serious problems for 
studying the heavy systems where the number of different configurations increases tremendously and almost always 
the algorithm remained stuck in a local minima. To deal with the more interesting large systems a sophisticated 
algorithm is needed which can handle the huge number of different configurations and finds the configuration with 
maximal binding energy in a reasonable amount of computational time. In addition, it should be able to overcome 
any type of local minima. 

We here present the details and technical aspects of such a new algorithm which is based on the simulated annealing 
method and is quite general in nature. The question addressed here requires an numerical approach which is not 
specific to the problem described. Apart from multifragmentation the energy minimization is needed, for example, 
in the nuclear structure calculations, in cluster radioactivity, in hadron physics etc. In nuclear cluster radioactivity, 
one is interested in the relative yields of different fragments which are emitted by the decaying nucleus. There, one 
assumes the isotopic distribution and the energy minimization is needed to find the most bound isobaric distribution 

Naturally, before one can talk about cluster formation, one first needs the phase space coordinates of the particles. 
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We here use the Quantum Molecular Dynamics (QMD) model || (as an event generator) to generate the time evolu- 
tion of the phase space coordinates of nucleons in a nuclear reaction. 



Our paper is organized as follow: The section II deals with the short description of the QMD model and a detailed 
description of the algorithm . The numerical realization of the algorithm is presented in section III and we summarize 
the results in section IV. 

II. The Formalism 

We here summarize shortly the Quantum Molecular Dynamics [QMD] model and then give the details of our new 
algorithm designed for multifragmentation. For more details on the QMD approach, we refer the reader to 13] . 

(i) The QMD approach: 

The QMD model is based on molecular dynamics and hence is an n-body theory which simulates the heavy ion 
reactions between 30 A ■ MeV to 1A ■ GeV on an event by event basis Jl0|]. Here each nucleus is represented by a 
coherent state of the form ( % =1 ) 

<t> a ( Xl ,t) = (J^J ' e -(*i-«<*-%r) 2 / 2i e <ft.(*i-*«) e -#. (1) 

The wave function has two time dependent parameter x a , p a . We fix the Gaussian width (L) to 1.08 fm 3 . In 
QMD calculations, nucleon a moves on a quasi-classical trajectory as obtained by a variation solution of the n-body 
Schrodinger equation: 

2ot = — + Vp a y^<V a p(x a ,xp,p a ,pp)>, (2) 

V 

Pa = -Vx a ^ < Vap(Xg,Vp,Pa,Pp) > ■ (3) 
P 

Here p a and x a are the centriods of the Gaussian wave functions in momentum and coordinate space which represent 
the nucleon a. The potential has the form M 

<V a p(x a ,Xp) >= Jd 3 Xld 3 X 2 <4>a4'l3\V(xi,X2)\4'a4>p>- (4) 

In addition, the nucleons interact via stochastic elastic and inelastic NN collisions. In principle, our approach to find 
the fragments is independent of the algorithm which generates the phase space coordinates. Therefore, QMD may 



be replaced by any other model ( like simple molecular dynamics model ||, Boltzmann-Uhling-Uhlenbeck model etc. 
p[ ) which is able to generate the phase space coordinates of the particles. Due to its n-body nature, the QMD model 
is more appropriate to study the fragment formation in heavy ion collisions than one body models. 



(ii) A Survey of Heavy Ion Reaction: 

During the simulation of the reaction, we store the phase space coordinates of all nucleons at several time steps. 
As the QMD model simulates the time evolution of nucleons , the stored phase space distribution is that of nucleons 
only. Our basic assumption is that in nature that configuration is realized which gives the largest binding energy. 
Therefore, a method has to be adopted to group the nucleons in free nucleons and fragments. The nucleons within 
a fragment will be bound by some binding energy. In a very simple model, one could consider the nucleons being a 
part of the same fragment if their centriods are closer than some spatial distance r max . This model is called minimum 
spanning tree [MST] method ||?] . One generally takes 2 < r max < 4. By definition, this method cannot address 
the fragment distribution during the violent phase of the reaction where whole nuclear matter is compressed and is 
confined to few fermis. The MST method at this time will give one single large fragment. More disturbing, the frag- 
ments detected by the MST method can contain nucleons with very large relative momenta. These fragments will be 
unstable and will decay after a while by fissioning or by emitting nucleons. To improve the model, a cut in momentum 
space has been also suggested recently by one of us and collaborators |?J . This cut (which limits the maximal allowed 
relative momentum of two nucleons in the same fragment) is quite effective in central collisions where most of the 
fragments are created during a reaction, but has no effect on the fragment distribution in peripheral collisions where 
the fragments are produced due to the decay of the spectator matter. 

If one combines the cuts in momentum and in coordinate space to a binding energy cut, one sees that several 
groups of nucleons are indeed not fragments, but a group of unbound nucleons which are close in spatial space. One 
has to follow the reaction for a long time until this group of nucleons decays in light and heavy fragments which are 
well separated in the coordinate space and can be detected with the standard MST algorithm. The critical time is 
generally assumed to be about 300 fm/c. 
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To give the reader a more clear picture, we simulated the reaction Au-Au at 600 MeV/nucl. and at impact param- 
eters of b = 3 and 8 fm, respectively, and display some key quantities in fig. 1 . The solid and dotted lines represent 
the reaction at 8 and 3 fm, respectively. The first row shows the evolution of mean density and of the collisions as 
a function of time. As expected, a higher density and collision number can be seen in central collisions compared to 
peripheral collisions. One also notices that the high reaction rate terminates at about 40-60 fm/c. Afterwards, we 
observe only collisions of nucleons in the same fragment. The second row shows the evolution of spectator ( filled 
circle) and participant ( filled triangle) nucleons. A participant nucleon is defined as a nuclcon which has undergone 
at least one collision. One sees that in central collisions 99% of nucleons have experienced a collision until 40 fm/c. 
As a results the directed transverse flow saturates as early as 40 -60 fm/c. The fig. 1(e) displays the evolution of the 
size of the largest fragment A max detected by the normal MST method with r max = 4 fm. We see one big fragment 
(consisting of 394 nucleons) at the time when the density is high. After about 120 fm/c we are able to find the "stable" 
fragment, which still decreases in size due to evaporation. Is this a realistic identification of the largest fragment? To 
answer this question, we applied a binding energy cut on the fragments detected in MST method. We first analyze 
the fragments with MST method and then pass all the fragments (with mass > 3) through an energy filter which 
recognizes a fragment only if it has at least a binding energy of 4 MeV/nucl. Otherwise it considers the MST fragment 
as a set of free nucleons. This approach is labelled as MST*. In both central and peripheral collisions, the largest 
fragment detected in MST is not a bound fragment at intermediate times. One gets properly bound fragment after 
about 120 fm/c after emitting the nucleons what lowers of course the binding energy. One should keep in mind that 
in peripheral collisions, one has two big (spectator) fragments and a fireball at mid-rapidity region without fragment. 



(iii) Simulated Annealing Clusterization Algorithm [SAC A] : 

Our new approach can be summarized as follows. We assume that : 

fl~| The nucleons from target and projectile are grouped into fragments (of any size) and into free nucleons. 
[2~| Though the nucleons inside a fragment can interact with each other, they do not interact with the nucleons 
from other fragments or free nucleons. 



3 That pattern of nucleons and fragments is realized in nature which gives the highest binding energy. 
To avoid that at intermediate times too many fragments are assumed (which finally break apart), we employ in 
addition a binding energy check. In order to form a fragment, the considered group of nucleons has to have a minimal 



binding energy given by 



(fa - Pn 1 }) 2 + m l - m <* + o V <*P( X <x> a : 1 



< L be x AT/, (5) 



with L(, e = -4.0 MeV if N* > 3 and Lt> e = otherwise. In this equation, N* is the number of nucleons in a fragment, 
Pprf is the center-of-mass momentum of the fragment. The binding energy criteria will make sure that no loosely 
bound fragments are formed in our approach. In reality these loosely bound fragments are not stable and decay during 
the reaction. The problem is that we have to find the most bound configuration among a huge number of possible 
patterns ( composed of nucleons and fragments). In order to cope with this complicated problem, we employ the 
simulated annealing technique and hence this algorithm is dubbed as " Simulated Annealing Clusterization Algorithm 
(SACA)" . 

One is tempted to start the search for the most bound cluster configuration by an iterative minimization method 
(also known as neighborhood search or local search). In this method, starting from a given configuration a new one is 
constructed. The new configuration is accepted only if it lowers the binding energy. The drawback of this procedure 
is that it may terminate at a local minimum. To improve this limitation, several modification can be imagined fl^] : 

1. To execute the algorithm for a large number of the initial configurations. This will finally allow to reach the 
global minimum. This is very time consuming. 

2. To use a algorithm which can jump over local minima and hence one can reach the global minima. This clearly 
depends strongly on the problem. Therefore its applications are limited. 

3: To generalize the iterative method so that the transitions which yields a higher binding energy are always 
accepted. In addition, the transitions which yield a lower binding energy are also accepted with a certain probability. 
This algorithm is known as simulated annealing method Jl2| . Its name is based on the fact that this algorithm is 
akin to the one used for cooling the solids. The simulated annealing method is a sequence of metropolis algorithms 
p3| with decreasing control parameter The control parameter § can be interpreted as a "temperature". For each 
Metropoliscity at a given temperature, we perform a sequence of steps until the binding energy does not change 
anymore. Each step is executed as follows: 

1: Given some initial configuration ip with energy f^,, a new configuration <p with energy is generated in the 
neighborhood of tp using a Monte-Carlo procedure. 

2: Let the energy difference between tp and (p is A£ = -f^,. 
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3.: If A£ is negative, the new configuration is always accepted. If A£ is positive, it is accepted with a probability 
exp(— A£/$). At the start, the control parameter •& is taken to be large enough for that most all attempted transitions 
are accepted. This is to overcome any kind of the local minima. After the binding energy remains constant, a gradual 
decrease in the control parameter z9 is made and the Metropolis algorithm is repeated. 

Note that there is no change in the coordinates of the nucleons. One should also note that employing this method 
does guaranty in the limit of infinite steps to reach the global minimum. Evidence that one reaches the ground state 
can be provided by obtaining the same fragment pattern for different starting configurations. 



To start with, a random configuration tp ( which consist of fragments and free nucleons) is chosen. The total energy 

associated with -0 configuration is given by 

n{ ( _ i N[ 

fc = EW ^ a ~ ^V 2 + m2 ~ ma + 2 S v <*p( x °" x p) 

Nl ( Nf 

a=l [ /3^q 

+ E \\l - pc N ^ 2 + m « - m « + \ E v °f>(**>*f>) \ 

a=l [ V " /3#a J f 

N L ( ^ N l 

+ ■ ■ ■ E \ y/&» ~ p ffl 2 + < - w « + 2 E Vopfroxp) 

a=l I fj^a 



Here is the number of nucleons in a fragment fi, P™} is the center of mass momentum of the fragment \x and 
V a p(x ai Xj3) is the interaction energy between nucleons a and (3 in a given fragment fi. Note that the total energy is 
the sum of the energies of individual fragments in their respective center of mass system. Therefore, ^ differs from 
the (conserved) total energy of the system because (i) the kinetic energies of fragments calculated in their center of 
masses and (ii) the interactions between fragments/free nucleons are neglected. At present a simple static interaction 
is implemented, but one can use the algorithm for arbitrary interactions. 



A new configuration is generated using Monte-Carlo procedure by either a) transferring a nucleon from some 
randomly chosen fragment to another fragment or by b) setting a nucleon of a fragment free or c) absorbing a free 
nucleon into a fragment. Let the new configuration ip be generated by transferring a nucleon from fragment v to 
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fragment [i. Then the energy of new configuration ip is given by: 



N( ( N( 

^ = E | y ®« - ^7) 2 + < - to « + 2 E ^(^. ^) 

a=l [ /3^a 

<+i r , <+l 

+ E | ^ ~ ^7+l )2 + ™« ~ ^ + 2 E Vap(x a ,X P ) 

+ -E< v^ a _ ^T )2 + < - m « + 2 E ^(^^f) 

a=l I /3#a 



Note that in this procedure, the individual energies of all fragments except for the donar fragment (v) and the receptor 
fragment {p) remain the same. The change in the energy from ip — ► <p is given by 

ac = C v -Ca- (6) 

Between the Metropolis algorithms, we cool the system by decreasing the control parameter A decrease in the 
temperature means that we narrow the energy difference which is accepted in a metropolis step. After many Metropo- 
lis steps, we should arrive at a minimum i.e. the most bound configuration. The problem is, however , that we usually 
arrive at a local minimum only. Between the local minimum, we find huge maxima. Let us give an example: Assume 
we have two fragments , but the most bound configuration would be one single fragment which combines both. Now 
each exchange of a single nucleons raises the binding energy and only the exchange of all nucleons at the same time 
lowers the total binding energy. This effect is well known in chemistry, where it is called Activation energy. In order 
to avoid this, we add , therefore, a second simulated annealing algorithm in which not anymore the nucleons are 
considered as the particles which are exchanged in each Metropolis step (like in the first simulated annealing), but the 
entities ( fragments or nucleons) obtained after the first step. This second stage of minimization is called fragment 
exchange procedure. This fragment exchange procedure is capable of overcoming any local minima. 



Note that even in this second stage of the minimization, the free nucleons can be exchanged as before. The total 



energy associated with any configuration W during second stage of iterations is given by 




a=l 



1 N Sl 

iPa - ) 2 + ml - m a + - Y V a p(x a , Xfj) 



-1/1 




Ns„ 



(fa - ) 2 + m 2 a - m a + - Y V a p(x a ,xp) 

P^a 

5 ! w ^ 

y (Pa - ^™ ) 2 + m l - ma + 2 X! ^(^a, a;/3) 

j W S „ 

(Pa - ^AT™ ) 2 + m a ~ m a + ^ E ^"/3 ( X " ' ^Z?) 



Here = X^=i* is the number of nucleons in a super-fragment S^. iVg is the number of nucleons in the i-th 
fragment contained in the super-fragment and Ng is the number of prc-fragments contained in the super-fragment 
Sfj,. The P^r™ is the center of mass momentum of the super fragment 5 M and V a p(x a , xp) is the interaction energy 
between nucleons a and [3 in a given super- fragment. Note that now the particle a interacts with its fellow nucleons 
in the same pre-fragment and also with the nucleons of other pre-fragments which are contained in a new given super 
fragment S^. 



Now the new configuration is generated using Monte-Carlo procedure by either a) transferring a pre-fragment from 
some randomly chosen super-fragment to another super-fragment or by b) setting a pre-fragment free or c) absorbing 
a single isolated pre-fragment into a super- fragment. Let us suppose that a new configuration $ is generated by 
transferring a pre-fragment i ( with mass N s ) from super-fragment v to super-fragment fx. The associated energy of 
new configuration $ reads as : 



'N Sl 

E 



N Sl 



~ 1 

(Pa - Ppfs 1 ) 2 + m «-™o + 7^ V a p(x a ,X ) 



+ ••■ 



+ 



E 

a=l 

E 

a=l 



P^a 



{Pa - P C £^ _ N i ) 2 + m 2 a - m a + | ^ V a0 {x a , xp) 



p^a 



(Po 



) 2 + ml - m a + - Y V af} (x a ,xp) 



Af s „ 



(Pa - ^AT™ n ) 2 + ™« - m Q + - Y Vap(x a , Xp) 



p^a 



The only difference between the particle and the fragment exchange procedure occurs for the bound nucleons. Now 
the bound nucleons cannot change their identity neither by being absorbed nor by becoming free. They will remain 
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bound in a pre- fragment. The pre- fragment itself can change its identity by either getting transferred to a new 
super- fragment, or be set free. As in the first stage, we calculate the energy difference between the new and the old 
configurations A£ and the metropolis procedure is continued till the most favored configuration is obtained. 

In summary, the actual procedure is as follows: We first simulate the nucleus-nucleus collision using the QMD 
model and store the phase-space coordinates of all nucleons at several time steps. At each stored time step, we apply 
the SACA to find the most bound configuration which consists of nucleons and fragments (of any size) . For a faster 
convergence of the algorithm, any cluster decomposition irrespective whether it fulfills the binding energy check (eq. 
H) or not is considered. Therefore, it is likely that several clusters may fail to fulfill eq. (||). At the end of the algorithm 
when the most bound configuration is found, we check the binding energy ( eq. 5) of each superfragment explicitly 
and mark all super-fragments violating this condition. The nucleons belonging to an 'inhibited' (marked) cluster are 
further on treated as free nucleons. The minimizing procedure of the simulated annealing mechanism is invoked again 
until a configuration is found where all fragments fulfill eq. (^). The heavy fragments are usually more bound than 
the lighter ones. We have carried out a detailed analysis and found that these are always the light fragments ( with 
masses 3 or 4 ) which at the end of the iterations are unbound or loosely bound. 

In the following, we discuss the numerical realization of the algorithm and present a detailed analysis of the influence 
of different parameters used in simulated annealing method. 

III. Numerical Realization: 

The simulated annealing algorithm has several parameters to be determined : the initial and the final value of the 
control parameter #, the number of metropolis steps to be executed at a given value of control parameter (i. e. length 
of Markov chain) , the decrease of the control parameter and the termination of the algorithm. This set of parameters 
is also referred as cooling schedule in the literature |l2| . One needs to choose the following parameters explicitly: 

1. : The initial value of the control parameter This will be referred as temperature. 

2. : The final value of the control parameter i}f [ i.e. the termination procedure]. 

3. : The length of the Markov chain M c h- 
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4.: A rule to fix the decrement in the control parameter a. 

Following , we use a so called simple cooling scheme and present the analysis of our extensive tests made for the 
collisions of Au-Au at 600 MeV/nucl. and at an impact parameter of 8 fm. We have also analyzed the results for the 
collisions of Pb-Pb (central) and Nb-Nb (central and peripheral) . The results of our analysis are independent of the 
masses of the colliding nuclei and also of the impact parameter. For our analysis we chose a conservative set of the 
above parameters and then try to find an optimized set of the parameters which yields the shortest computational 
time. We use the following set of parameters if not stated otherwise: 

The initial temperature is taken to be 4 MeV. The length of Markov chain is taken to be 70??. ;rj being the 
number of entities at the start of the minimization. After each markov chain [ = 70?/] , the temperature is decreased 
using a simple law : 

with a — 0.95. Finally, the algorithm is terminated if there is no change in the binding energy for a large number of 
iterations [ = 6O77] . The details of each of these parameters will be presented the following paragraphs. 

[Tj The Initial Configuration: We have to choose a random initial distribution initially to evoke the simulated an- 
nealing minimization. In our procedure, we distribute the nucleons (of the two colliding nuclei) into few cells. The 
transfer of nucleons is allowed among these cells. Naturally, the final outcome should be independent of the number 
of cells we choose. In fig. 2, we present the outcome of a single QMD event when exposed to SAC A with a different 
number of initial cells. The displayed reaction is of Au-Au at 600 MeV/nucl. and impact parameter of 8 fm. Here 
we vary the number of cells between two and 394 ( that is by treating each nucleon as a free particle). We see that 
the variation in the cell number does not affect the final fragment distribution. At zero fm/c, the simulated annealing 
method finds two nuclei ( i.e. the projectile and target) which shows the validity of the annealing method. One also 
notices that the binding energy of the system remains constant between 40 fm/c and 200 fm/c. In other words, the 
most bound configuration found in SACA at 40 fm/c and 200 fm/c is approximately the same. 

In fig. 3, we display the evolution of most bound configuration using the two extremes : 2 cells and all particles 
free at 0, 40, 120 and 200 fm/c, respectively. Note that the high density phase is reached around 40 fm/c. Between 
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120 and 200 fm/c, one should not expect much change as the reaction is already finished and there is only some 
rearrangement of the nucleons in the fragments. In case of 2 cells ( two initial clusters) the algorithm first breaks each 
of the clusters into large number of free nucleons (because free nucleons have zero energy) and some light fragments. 
After several hundred thousands iterations, it starts rearranging the nucleons into bound fragments. It is interesting 
to note that after some initial differences, the evolution of most bound configuration is quite the same in both cases. 



As stated in the algorithm section, we choose the new configuration (ip or $) by transferring a nucleon/pre- fragment 
from one fragment/superfragment to another. These fragments are chosen by Monte-Carlo method. It would be in- 
teresting to study the effect of different Monte-Carlo procedures (applied in SACA) on the fragment distribution. 
The different Monte-Carlo procedures can be generated using different random numbers. In figure 4, we display 
the fragment distribution ( i.e. the largest fragment A rnax , the number of free nucleons and of intermediate mass 
fragments IMF's A > 5 ) and the energy of the most bound configuration at three different times i.e. at zero fm/c, 40 
fm/c and 200 fm/c, respectively, for the different iterations. We see, as it should be, a almost complete independence. 

|~ji~| . The Initial Value di : One of the key features of the Metropolis algorithm is that it also accepts the transitions 
which increases the cost function i.e. the energy. Therefore, the initial value di should be such that the most of the 
attempted transitions are accepted during first iterations. In other words, exp(— A£/?^) ~ 1. A practical way to 
implement the sequence of ?Vs is given by Johnson et.al. J3. There the average increase in the energy over large 
number of iterations AC is related with $i by 

acceptance ratio \ = exp {— A£/^i} . (7) 

i.e. 

Generally, the acceptance ratio \ should be close to 1. The choice of the value di (or the temperature) depends very 
strongly on the problem at hand. It should be kept in mind that a very large value of di will lead to huge computa- 
tional time whereas a very small value will lead to less attempts which are accepted by the Metropolis algorithm and 
consequently which may lead to a wrong final distribution. In fig. 5, we show the same reaction as reported in fig. 
1-4, but at different values of initial temperature di. Here the other parameters are kept unchanged. The variation 
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of j?j between 1 MeV and 500 MeV has no effect on the fragment distribution at zero fm/c. We have just two gold 
nuclei initially. On contrary, one can see some differences at 40 fm/c. A very small value of i?j (< 3-4 MeV) leads to 
a heavier A max (95) compared to the average A max ( w 42) and as a result fewer IMF's and nucleons are emitted. 
Similar conclusions can be drawn at 200 fm/c. A very low value of i?j apparently freezes the initial configuration. The 
results are more stable for > 3MeV. Therefore, we choose the #j= 5 MeV. 



In fig. 6, we display the evolution of the most bound configuration in Au-Au reactions using di = 5 MeV and 
500 MeV. We see that when one iterates the reaction with very large initial temperature (=500 MeV), almost all 
attempted transitions are accepted in the Metropolis algorithm. With a moderate value of the temperature i?j = 5 
MeV, only some selected configurations are accepted. The minimization of the energy with #j= 500 MeV results in 
the vibration around the same fragments for very long time. An (unnecessary) large value of i?j does not help to 
establish an early equilibrium. In contrary, one needs huge computation time (in terms of iterations) to find the most 
bound configuration. The same can be achieved with moderate value of #j = 5 MeV with far less costs. 



The Length of the Markov Chain: M c h : Here we fix the control parameter $ and execute the algorithm for 

a fixed number of Metropolis steps. We construct a sequence of fragment configurations Q = {?/>, <p, , $}. One 

should note that here we have an initial configuration ip and new configuration (p(= ip + 1) is generated by a random 
matrix. Thus, the number of iterations, and hence the length of Markov chain should be long enough to ensure an 
equilibrium. In fig. 7, we show the results with Markov chains of different length. The length of the markov chain n is 
displayed in the units of the total number of the nucleons ( prefragments) present at the beginning of the minimization. 
After about 40 r\ the results are quite stable. For a smaller values of 77, there are fluctuations in the results and in 
addition, SACA overestimates the size of A max and underestimates consequently the IMF's production. We fix the 
length of the markov chain rj~ 40. The effect of different M c ^'s on the evolution of the most bound configuration is 
shown in fig. 8 where the evolution of the fragment's multiplicity is plotted as a function of the iterations for two 
values of Markov's chains i.e for M c h = 40 77 and 450 77, respectively. Note that the number of iterations is = rj times 
the number of temperature steps. We see that the initial evolution is quite the same in both cases, but a smaller value 
of M c h needs less iterations than a longer one to arrive at same final value. This is easy to understand. The main aim 
of iterating over large number of iterations with same 1? is to establish the quasi-equilibrium. Once an equilibrium is 
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established, there will be no further improvements in the cost at same therefore, a smaller value of M c h leads to 
same result as that with largest M c h- 

| iv | . Decrement in the Control Parameter : The minimization is started with a relative large temperature $j . 
Then the temperature is decreased in steps after a quasi-equilibrium is established for each temperature. Apparently, 
a larger decrement in the temperature will lead the defect to be frozen i.e. any configuration which may or may not 
be the most bound can freeze whereas a very small decrement will need huge computational time. The decrement 
should be in such a way that the length of the Markov chain M c h is as small as possible and thus after a new change 
of control parameter i?, the quasi-equilibrium should be re-established as soon as possible. We here follow the simple 
rule for the decrement factor a. 

= o~ ' $i (9) 

The value of a varies in the literature between 0.5 to 0.95 jl2), fll4|| , fll(| . The effect of different decrement factors 
a is displayed in fig. 9. Here all other parameters are kept the same as discussed at the beginning. One can see that 
a very small value of a overestimates the size of A max and underestimates the IMF production. We fix the value of a 
to 0.85. The comparison of two simulations resulting from the decrement factor er= 0.85 and 0.98 is displayed in fig. 
10. Here we see that the two different values gives the same cooling result but for the larger value of a many more 
iterations are necessary. 



The Final Value The termination procedure used in the literature varies from problem to problem and 
also from author to author. We fix the termination by two different controls. 

1. Either we stop the calculations if the control parameter $ has reached a very small value where no further 
transition can be expected. For the present calculation, we take i)f — 10~ 10 MeV. 

2. Or we terminate the algorithm if there is no change in the configuration over a large number of attempted 
iterations.. Following the rule used to fix the length of the Markov chain, we choose the length for termination Iterm 
in terms of r\ which represent the number of iterations in the units of M c h- In fig. 11, we display the effect of the 
variation in Iterm on the fragment distribution. We find that different termination lengths do not affect the results. 
The effect of Iterm on the evolution of the most bound configurations is shown in fig. 12 where the evolution as a 
function of iterations is displayed for two values of Iterm- i.e for Iterm = 5 r\ and 120 77, respectively. The different 
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termination values have a very small effect on the fragment structure. Therefore, we fix l te rm to 35 rj. 
In above paragraphs, we have discussed in detail the influence of different choices of the parameters which determine 
the simulated annealing method. One should note that once these parameters are chosen, the simulated annealing 
method is completely determined and it is a complete self-iterative method. 

In the further discussion, the set with conservative parameters ( i.e. with di — 500 MeV, M c h — 450?7, a — 0.98 and 
kerm = 120?y) is called as Si whereas the set with the most economical parameters ( i.e. with di = 5 MeV, M c h = 40??, 
er = 0.85 and kerm — 35 -q) is called as S ec . 

The crucial test of the algorithm is its application to a single nucleus in its ground state. In principle one should get 
a single nucleus at the end. But the results can be different in reality To see the importance of nucleon and fragment 
exchange procedured in SACA, we first turned off the fragment exchange part of the algorithm ( i.e. 2nd stage of 
the algorithm). As expected, the transfer of single nucleons terminates in local minima and as a result, one finds 
several fragments in the ground state of a nucleus. Naturally, the global minima is a single nucleus. The energy gain 
in a single nucleon transfer is not enough to overcome the huge energy barrier. This energy barrier can be overcome 
only by allowing the collective transfer of the nucleons i.e. of fragment as such. When we turned on the fragment 
exchange procedure of the algorithm, we could overcome the local minima and we find single nucleus as most bound 
configuration. In fig. 13, we show the evolution of the fragments as a function of the iterations for different single 
nuclei 20 Ne, 40 Ca, 93 Nb and 208 Pb, respectively using S ec . Starting points are the (ground) state nuclei as generated 
by the QMD. In all cases, the SACA finds the single nucleus at the end of the iterations as it should. One also notices 
that the lighter nuclei need less iterations to find the most bound configuration. We find « 8,000,12,000,62,000 
and 280, 000 iterations are necessary to find the ground state for 20 A^e, 40 Ca, 93 Nb and 208 Pb, respectively. One also 
notices that the increase in the number of necessary iterations is not a linear function of the masses of nuclei. The 
energy of the configurations is displayed in fig. 14. We notice that one has a positive energy at the beginning which 
is decreased by breaking the cells into large number of nucleons/fragments. After a large number of iterations, one 
finally reaches the most bound configuration. 

In fig. 15, we display the multiplicity (averaged over 20 events) of different fragments obtained using Si and S ce , 
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respectively. Here Au-Au collision is carried out at impact parameter of 8 fm. We see that both sets of parameters 
give a similar evolution of the reaction. One should note that the minimization with S ec needs much less computing 
time as compared to Si. Our algorithm is able to detect the fragment distribution as early as 50-60 fm/c. From fig. 
1, one notices that the density is maximum at this time. This very early identification of fragments in SACA is very 
promising because it means that the fragments may give insight into hot and dense nuclear matter. 

The annealing algorithm can be made faster if some pre-information is feed into the algorithm. Naturally, the 
nucleons which are very far away in spatial or in momentum space will not lower the energy if one combines them 
as fragment. We applied a cut in spatial and in momentum space to sort out those distant nucleons. We took a 
minimal spatial distance between two nucleons of 10 fm and a relative momentum of 200 McV/c. In other words, we 
first break the whole system into fragments using these conditions and each of these fragments are then subjected to 
SACA. We found that the results are the same as before but the algorithm is about 10 times more faster. 



IV. Summary: 



Summarizing, based on the simulated annealing method, we have presented the details of a new algorithm devel- 
oped to study multifragmentation of heavy ion collisions. We have carried out an extensive survey of the different 
parameters which are crucial for the success of the method. Based on our calculations, a set of parameters is suggested 
for the algorithm which makes the algorithm very fast and accurate. This algorithm can detect the fragments as early 
as 40-60 fm/c i.e. at time when density is relative high and the interactions between fragments are still going on. It 
can not only give insight into the hot and dense nuclear matter, but at the same time makes it possible to apply the 



full in-medium G-matrix approach |17| to study the multifragmentation which was not possible due to the emission 
of nucleons after 70-80 fm/c [Q. A brief outcome of the results was presented in ref. |lj| and a detailed physical 
interpretation of our results for various reactions will be presented elsewhere p0| . The algorithm is very general and 
may serve for every problem in which the most bound configuration has to be found. 
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Figure Captions: 



Fig. 1. Evolution of Au-Au collisions at incident energy of 600 MeV/ nucl. using a soft equation of state. The 
results at b = 3 and 8 fm are displayed, respectively, by dotted and solid lines. Fig. 1(a) displays the mean density 
whereas the rate of collision is shown in Fig. 1(b). The evolution of the spectators ( filled circle) and the participants 
(filled triangle) is are presented in Fig. 1(c). Fig. 1(d) shows the time evolution of the transverse flow of the nucleons. 
Here, we do not consider the formation of fragments. The evolution of the largest mass A max formed within MST 
and MST with binding energy check (MST*) are displayed in fig. 1(c) and Fig. 1(f), respectively. 

Fig. 2. The heaviest fragmncnt A max , the emitted nucleons, the multiplicity of fragments with mass A > 5 and 
the total energy associated with the configuration is displayed as a function of the cell number. Here the results at 
, 40 and 200 fm/c are represented, respectively, by filled circle, open square and filled triangle. The displayed results 
are for a single event generated using QMD model. 

Fig. 3. The evolution of the most bound configuration as a function of the iterations. Here we display the results 
at four times i.e. at 0, 40, 120 and 200 fm/c, respectively. The solid and dotted lines represents the results obtained 
with cells = 2 and and 394, respectively. 

Fig. 4 Same as fig. 2, but as a function of Monte -Carlo procedures. 

Fig. 5 Same as fig. 4, but as a function of the initial temperature Here the arrow shows the value of the 
parameter chosen for the optimized set of parameters. 

Fig. 6 Same as fig. 3, but with temperature #j= 5 MeV ( solid line) and 500 MeV ( dotted line), respectively. 

Fig. 7 Same as fig. 5, but as a function of the length of the Markov chain M c h- 

Fig. 8 Same as fig. 3, but with M c h = 407/ and 450r7, respectively. 
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Fig. 9 Same as fig. 5, but as a function of the decrement factor a. 
Fig. 10 Same as fig. 3, but with a — 0.85 and 0.98, respectively. 
Fig. 11 Same as fig. 5, but as a function of termination length Iterm- 
Fig. 12 Same as fig. 3, but with l te rm = 35r; andl20?7. 

Fig. 13 Same as fig. 3, but the evolution of single nuclei Ne, Ca, Nb and Pb, respectively. 
Fig. 14 Same as fig. 13, but the energy of the system as a function of the iterations. 

Fig. 15 The time evolution of the Au-Au at 600 MeV/nucl and at an impact parameter of 8 fm using a soft 
equation of the state. Here we display the results which arc averaged over 20 events. The results obtained using Si 
and S ec are shown , respectively, by the filled circles and open squares, respectively. 
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